#include "CHBDYPLINECalculator.h"
#include "Utility/Algebra/Tools.h"
#include "DataStructure/Input/Solver/ThermodynamicsConstant.h"

#include "DBManager/IService/Input/Node/INodeService.h"
#include "DBManager/IService/Input/Element/IElementService.h"
#include "DBManager/IService/Input/Subcase/ISubcaseService.h"
#include "DBManager/IService/Input/Material/IThermalMaterialService.h"
#include "DBManager/IService/Input/Property/IPCONVService.h"
#include "DBManager/IService/Input/Load/ICONVService.h"
#include "DBManager/IService/Input/Constraint/ITEMPService.h"
#include "DBManager/IService/Input/Constraint/ISPCService.h"
#include "DBManager/IService/Input/Property/IPHBDYService.h"
#include "DBManager/IService/Input/Load/IQBDY3Service.h"

#include "SGFem/Pre/Assembly/Assembly.h"
#include "SGFem/Pre/Assembly/ModelInfoQuery.h"


using SG::Algebra::Point;
using SG::Algebra::Vector3D;
using SG::Algebra::Matrix;
using SG::Algebra::length;
using SG::DataStructure::Common::Status;
using SG::DataStructure::Common::Real;
using SG::DataStructure::Common::TOLERANCE;
using SG::DataStructure::Common::Id;
using SG::DataStructure::FEM::ThermalMaterial;
using namespace SG::DataStructure::FEM;
using namespace SG::DBManager;
using SG::FEMSolver::AssemblyEigen;

using SG::FEMSolver::Element::CHBDYPLINECalculator;
std::size_t CHBDYPLINECalculator::s_nNode = 2;
std::size_t CHBDYPLINECalculator::s_nDim = 2;


void CHBDYPLINECalculator::Initialize (SG::DataStructure::FEM::ElementTypes type)
{
    m_type = type; 
}

void CHBDYPLINECalculator::ComputeConvStiffness
    (
    const std::vector<Point>& points, 
    const SG::DataStructure::FEM::ThermalMaterial& material,
    const Real& areafactor,
    _OUT SG::Algebra::Matrix& ke2
    )
{
    Real h =  material.m_H;
    Real lenij = length(points[0], points[1]);
    ke2 = Matrix (2, 2, {2.0, 1.0, 1.0, 2.0});
    ke2 = ke2 * (areafactor * h  *lenij / 6.0);
}

void CHBDYPLINECalculator::ComputeConvStiffness
(
    const std::vector<Point>& points, 
    const std::vector<Real>& h,
    const Real& areafactor,
    _OUT SG::Algebra::Matrix& ke2
)
{
    Real hsum = (h[0] + h[1]) / 2.0;
    Real lenij = length(points[0], points[1]);
    ke2 = Matrix (2, 2, {2.0, 1.0, 1.0, 2.0});
    ke2 = ke2 * (areafactor * hsum *lenij / 6.0);
}

void CHBDYPLINECalculator::ComputeRadStiffness
    (
    const std::vector<Point>& points, 
    const SG::DataStructure::Common::Real& Emissivity,
    const std::vector<SG::DataStructure::Common::Real>& T,
    const SG::DataStructure::Common::Real& Text,
    const Real& areafactor,
    _OUT SG::Algebra::Matrix& ke4
    )
{
    if (T.size() != 2) {
        // 丢出错误并终止程序
        throw std::runtime_error("Error: The radiating boundary temperature points of the line element should be 2.");
    }
    
    std::vector<Real> hr(s_nNode);
    for (std::size_t i{0}; i < T.size(); ++i) 
    {
        hr[i] = STEFAN_BOLTZMANN_NUMBER * Emissivity * (T[i] * T[i] + Text * Text) * (T[i] + Text);
    }

    ComputeConvStiffness(points, hr, areafactor, _OUT ke4);
}


void CHBDYPLINECalculator::ComputeConvLoad
    (
    const std::vector<Point>& points, 
    const SG::DataStructure::FEM::ThermalMaterial& material,
    const SG::DataStructure::Common::Real& Text,
    const Real& areafactor,
    _OUT SG::Algebra::Matrix& p3
    )
{
    Real h =  material.m_H;
    Real lenij = length(points[0], points[1]);

    p3 = Matrix (2, 1, {1.0, 1.0});
    p3 = p3 * (areafactor * h * Text * lenij / 2.0);
}

void CHBDYPLINECalculator::ComputeConvLoad
    (
    const std::vector<Point>& points, 
    const std::vector<Real>& h,
    const SG::DataStructure::Common::Real& Text,
    const Real& areafactor,
    _OUT SG::Algebra::Matrix& p3
    )
{
    Real hsum =  (h[0] + h[1]) / 2.0;
    Real lenij = length(points[0], points[1]);

    p3 = Matrix (2, 1, {1.0, 1.0});
    p3 = p3 * (areafactor * hsum * Text * lenij / 2.0);
}

void CHBDYPLINECalculator::ComputeHeatFluxLoad
    (
    const std::vector<Point>& points, 
    const SG::DataStructure::Common::Real& q0,
    const Real& areafactor,
    _OUT SG::Algebra::Matrix& p2
    )
{
    Real lenij = length(points[0], points[1]);

    p2 = Matrix (2, 1, {1.0, 1.0});
    p2 = p2 * (areafactor * q0 * lenij / 2.0);
}

void CHBDYPLINECalculator::ComputeHeatFluxMatries 
    (
        const SG::DBManager::DBServiceFactorySP dbServiceFactory, 
        const SG::DataStructure::FEM::CHBDYPData eleData,
        const SG::DataStructure::Common::Real& q0,
        const SG::DataStructure::Common::Id sid
    )
{
    auto& assembly = AssemblyEigen::GetInstance ();
    
    static const DofSet DOFs {DOFDirections::DOF_Ux};
 
    PHBDYData phbdyData;
    std::vector<std::size_t> nodeIndexes(s_nNode);
    std::vector<Id> dispCoordIds(s_nNode);
    std::vector<int> dofGlobalIndex(s_nNode);
    std::vector<Point> points(s_nNode);
  
    const auto& nodeService = dbServiceFactory ->get<INodeService>();      
        
    //单元特性数据
    dbServiceFactory ->get<IPHBDYService>()->FindById(eleData.m_pId, _OUT phbdyData);

    // Node coordinate
    ModelInfoQuery::QueryNode (nodeService, eleData.m_g, s_nNode, _OUT points, _OUT nodeIndexes, _OUT dispCoordIds);

    assembly.GetGlobalDOFId (nodeIndexes, DOFs, _OUT dofGlobalIndex);
    
    //热流荷载贡献
    Matrix p2(s_nNode, 1);
    ComputeHeatFluxLoad(points, q0, phbdyData.m_AF, _OUT p2);
    assembly.AccumulateLoad(sid, dofGlobalIndex, p2);
    
}

void CHBDYPLINECalculator::ComputeConvMatries 
    (
        const SG::DBManager::DBServiceFactorySP dbServiceFactory,
        const SG::DataStructure::Common::Id sid
    )
{
    auto& assembly = AssemblyEigen::GetInstance ();
    
    static const DofSet DOFs {DOFDirections::DOF_Ux};

    CHBDYPData eleData;   
    PHBDYData phbdyData;
    PCONVData property;
    CONVData convData;
    
    //TEMPData tempData;
    //SPCData spcData;
    //QBDY3Data heatFluxData;
    ThermalMaterial material;
    
    //std::vector<QBDY3Data> heatFluxDataArray;
    std::vector<TEMPData> tempDataArray;
    
    std::vector<std::size_t> nodeIndexes(s_nNode);
    std::vector<Id> dispCoordIds(s_nNode);
    std::vector<int> dofGlobalIndex(s_nNode);
    std::vector<Point> points(s_nNode);
    
    const auto& chbdypService = dbServiceFactory ->get<IElementService>();  
    const auto& nodeService = dbServiceFactory ->get<INodeService>();      
    const auto& convService = dbServiceFactory ->get<ICONVService>();
    const auto convSize = convService-> Size();
    for (std::size_t index{0}; index < convSize; ++index)
    {
        /** 从数据库获取对应的对流数据 */
        convService->FindByIndex (index, _OUT convData);
        
        if (Status::STATUS_Failure == chbdypService->FindById (convData.m_id, _OUT eleData))
        {
            continue;
        }
        
        //单元特性数据
        dbServiceFactory ->get<IPHBDYService>()->FindById(eleData.m_pId, _OUT phbdyData);

        // Node coordinate
        ModelInfoQuery::QueryNode (nodeService, eleData.m_g, s_nNode, _OUT points, _OUT nodeIndexes, _OUT dispCoordIds);

        assembly.GetGlobalDOFId (nodeIndexes, DOFs, _OUT dofGlobalIndex);
        
        //获取对流换热系数和热交换温度
        dbServiceFactory ->get<IPCONVService>()->FindById (convData.m_PCONID, _OUT property);
        dbServiceFactory ->get<IThermalMaterialService>()->FindById (property.m_MID, _OUT material);

        //对流刚度矩阵贡献
        Matrix ke2(s_nNode, s_nNode);
        ComputeConvStiffness (points, material, phbdyData.m_AF, _OUT ke2);
        assembly.AccumulateStiffness (dofGlobalIndex, ke2);

        //对流荷载贡献
        Matrix p3(s_nNode, 1);
        Real Text;
            
        dbServiceFactory ->get<ITEMPService>()->FindBySId(sid, _OUT tempDataArray);
        for (auto& tempData : tempDataArray)
        {
            if (tempData.m_G == convData.m_TA[0])
            {
                Text = tempData.m_T;
                break;
            }
        }
        ComputeConvLoad (points, material, Text, phbdyData.m_AF, _OUT p3);
        assembly.AccumulateBCLoad(sid, dofGlobalIndex, p3);
    }
}

void CHBDYPLINECalculator::ComputeElementMatries (const SG::DBManager::DBServiceFactorySP dbServiceFactory)
{
    auto& assembly = AssemblyEigen::GetInstance ();
    
    static const DofSet DOFs {DOFDirections::DOF_Ux};

    CHBDYPData eleData;   
    PHBDYData phbdyData;
    PCONVData property;
    CONVData convData;
    // TEMPData tempData;
    // SPCData spcData;
    QBDY3Data heatFluxData;
    ThermalMaterial material;
    // std::vector<CONVData> convDataArray;
    std::vector<QBDY3Data> heatFluxDataArray;
    std::vector<TEMPData> tempDataArray;
    
    std::vector<std::size_t> nodeIndexes(s_nNode);
    std::vector<Id> dispCoordIds(s_nNode);
    std::vector<int> dofGlobalIndex(s_nNode);
    std::vector<Point> points(s_nNode);
    
    const auto& chbdypService = dbServiceFactory ->get<IElementService>();  
    const auto& nodeService = dbServiceFactory ->get<INodeService>();      
    const auto eleSize = chbdypService->GetElementSize (m_type);
    for (std::size_t index{0}; index < eleSize; ++index)
    {
        /** 从数据库获取对应的单元数据 */
        chbdypService->FindByIndex (index, _OUT eleData);
        
        //单元特性数据
        dbServiceFactory ->get<IPHBDYService>()->FindById(eleData.m_pId, _OUT phbdyData);

        // Node coordinate
        ModelInfoQuery::QueryNode (nodeService, eleData.m_g, s_nNode, _OUT points, _OUT nodeIndexes, _OUT dispCoordIds);

        assembly.GetGlobalDOFId (nodeIndexes, DOFs, _OUT dofGlobalIndex);
        
        Id sid = 1; //这个好像是用spc来确定响应的id
        //获取对流换热系数和热交换温度
        if (Status::STATUS_Success == dbServiceFactory ->get<ICONVService>()->FindById(eleData.m_id, _OUT convData))
        {
            dbServiceFactory ->get<IPCONVService>()->FindById (convData.m_PCONID, _OUT property);
            dbServiceFactory ->get<IThermalMaterialService>()->FindById (property.m_MID, _OUT material);

            //对流刚度矩阵贡献
            Matrix ke2(s_nNode, s_nNode);
            ComputeConvStiffness (points, material, phbdyData.m_AF, _OUT ke2);
            assembly.AccumulateStiffness (dofGlobalIndex, ke2);

            //对流荷载贡献
            Matrix p3(s_nNode, 1);
            Real Text;
            
            dbServiceFactory ->get<ITEMPService>()->FindBySId(sid, _OUT tempDataArray);
            for (auto& tempData : tempDataArray)
            {
                if (tempData.m_G == convData.m_TA[0])
                {
                    Text = tempData.m_T;
                    break;
                }
            }
            ComputeConvLoad (points, material, Text, phbdyData.m_AF, _OUT p3);
            //assembly.AccumulateLoad(sid, dofGlobalIndex, p3);
            assembly.AccumulateBCLoad(sid, dofGlobalIndex, p3);
            //std::cout<<"convection test!"<<std::endl;
        }
        

        // sid = 2; //要把loadcase传进来
        // //获取热流参数
        // if (Status::STATUS_Success == dbServiceFactory ->get<IQBDY3Service>()->FindBySId (sid, _OUT heatFluxDataArray))
        // {   
        //     for (auto& heatFluxData : heatFluxDataArray)
        //     {
        //         if (eleData.m_id == heatFluxData.m_EID)
        //         {
        //             //热流荷载贡献
        //             Matrix p2(s_nNode, 1);
        //             Real q0 = heatFluxData.m_Q0;
        //             ComputeHeatFluxLoad(points, q0, phbdyData.m_AF, _OUT p2);
        //             assembly.AccumulateLoad(sid, dofGlobalIndex, p2);
        //         }
        //     }
        //     continue;
        // }
    }
}

// 后处理数据
void CHBDYPLINECalculator::DataRecovery (const SG::DBManager::DBServiceFactorySP dbServiceFactory, int subcaseId, SG::DataStructure::Common::Real time, const std::vector<SG::DataStructure::Common::Real>& globalDisplacement)
{
}

// 计算全局坐标系下单元应力刚度矩阵
void CHBDYPLINECalculator::ComputeStressStiffness (const SG::DBManager::DBServiceFactorySP dbServiceFactory, int staticSubId)
{
}
